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The nature of phase boundaries in the QCD phase diagram has not been satisfactorily explored 
by experiments. Based on the Ginzburg-Landau free energy with a spatially inhomogeneous term 
as a function of a scalar order parameter, it is possible to determine spatial correlation lengths 
from measured two point correlation functions in general, where a spatially dependent multi- 
plicity fluctuation from the mean value is taken as the order parameter. The divergence of the 
correlation lengths can be a robust signature of a critical system as well as the disappearance of 
< qq > condensations. Simultaneous observations of such observables can probe the nature of 
the phase boundary conclusively. In this letter, a present result from the PHENIX experiment 
will be reported. We will focus whether a critical behavior of the phase transition exists or not 
by searching for increases of pseudo-rapidity correlation lengths as a function of the number of 
participant nucleons with dominantly low p, charged particles where instant fluctuations such as 
high p t jets are expected to be negligible. The result shows a non monotonic increase in the 
correlation lengths which can be a symptom of the critical behavior. 
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1. Introduction 

Many observations at RHIC indicate the formation of the strongly interacting dense medium 
with partonic natures[[j]]. Nevertheless, information on the phase boundaries have not been quan- 
tified. Although lattice QCD calculations as well as model calculations predict the existence of 
a tricritical end-point, none of them reaches a quantitative agreement each other on the location 
of the end-point[Q]. Hence it is crucial to pin down the end-point in the QCD phase diagram and 
investigate the transition order by experiments. Before pinning down the end-point as an ultimate 
goal, it is important to establish ways to determine critical systems in general. As one of such 
observables, increases of spatial correlation lengths as a function of system energy density can be 
a robust signature to determine critical systems whatever the transition order is. 

In the following sections the relation between the phase transition and the density fluctuations 
is briefly reviewed and the experimentally accessible observables are advocated. Based on the 
observables, Au+Au collision events taken by the PHENIX detector^] at y/s^ = 200GeV have 
been analyzed and the present results are summarized by focusing on whether critical behaviors of 
the phase transition exist or not as a function of the number of participants N p which reflects the 
system energy density [Q]. 

2. Density fluctuation and phase transition 

In order to relate the density fluctuations with the phase transition in the simplest form, 
Ginzburg-Landau(GL)[|5|] theory with the Ornstein-Zernike picture^] for a scalar order parame- 
ter is briefly reviewed. The first attempt to apply the free energy discussion to nucleus-nucleus 
collisions can be found in [^]. GL describes the relation between a free energy density / and an 
order parameter <p as a function of system temperature T. By adding a spatially inhomogeneous 
term (V0) 2 and an external field h, the general form is described as follows; 

f(T,<j>,h) = /o (r) + iA(r)(V0) 2 + 

^a{T)$ 2 + h(f> 4 + h(j) (2.1) 

where terms with odd powers are neglected due to the symmetry of the order parameter and the 
sign of b is used to classify the transition orders; b < for first order, b > for second order and 
b = for the critical point. Since the order parameter should vanish above critical temperature T c , it 
is natural for the coefficient a{T) to be expressed as a(T) = ao(T — T c ), while b is usually assumed 
to be constant in the vicinity of T c . 

In the following analysis, the order parameter corresponds to the multiplicity density fluctu- 
ation from the mean density as a function of one dimensional rapidity point y, which is defined 
as 

Hy)=p(y)-(p) (2-2) 

where a pair of brackets is an operator to take the average. With the Fourier expansion of the 
density fluctuation (j) (y) = Y,k fye where k is wave number, one can express the deviation of the 
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free energy density AF /Y due to spatial fluctuations from the equilibrium value as 

AF/Y = ^J(f-f )dy=^\(t> k \ 2 (a(T)+A(T)k 2 ) (2.3) 

where Y is the total rapidity range corresponding to one dimensional volume and up to the sec- 
ond order terms are taken into account as an approximation in the vicinity of the critical point in 
Eq.([Q|). Given the free energy deviation, one can obtain the statistical weight w for fluctuation 
ty{y) to occur in a given temperature T 

w(<p(y))=Ne- AF/T . (2.4) 

Therefore the statistical average of the square of the density fluctuation with the wave number k is 
described as 



(W 2 } = J^M 2 w\Yhe ik >\d<$> k 



- NT 1 

~^a(T)+A(T)k 2 ' ( } 

Experimentally observable two point density correlation function can be related to the statis- 
tical average of the square of the density fluctuation. With a density p(yi) for a given sub volume 
dyt, the general two point density correlation G2 is expressed as 

<h(yi,yz) = ((pW- (P))(p(y 2 ) - <p))> (2.6) 

where the case that 1 coincides with 2 is excluded to simplify the following discussion. Multiplying 

e -iky = £-*%2-)>i) to t h e s j(j es Q f Eq^l^) and integrating over sub volume dy\ and dy% gives 

Y J G 2 (y)e- ik >'dy = (| J (p(y) - {p^ dy\ 2 ) = <|«^| 2 ). (2.7) 
From Eq.(|2~3|) and (|2.7|), G2 can be obtained by the inverse Fourier transformation of (|0/t| 2 )- There- 



fore in the one dimensional case G2 is described as 



G2(y) = ^^^(T)e- m(T \ (2-8) 
where a correlation length t, (T) is inuoduced, which is defined as 

In general, a singular behavior of t, (T) as a function of T indicates the critical point of the phase 
transition. 



The wave number dependent susceptibility can also be defined from Eq.(|2.1[) and ( |2.3| ) as 
follows 



_/aV\ (dhy l _( d 2 (AF/Y) \ 



a Q {T-Tc)(\+k^(TY 



1 . (2.10) 
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In the case of the static limit of k = 0, the susceptibility can be expressed as 

1 2Y 2 

*-° = M^=AT 5(T)G2(0) - (2 - H) 



From Eq.( |2.9[ ) and ( 2.1 1[ ), simultaneous observations of singularities or discontinuous behaviors in 



both E,{T) and Xk=oT x %(T)G2(0) at the same temperature indicate that the phase transition is 
consistent with the second order transition by assuming continuous variations of T. 

3. Experimental observables 

In the following analysis the density fluctuation is discussed via the charged particle multiplic- 
ity distributions as a function of the pseudo-rapidity interval of 5t] for each collision centrality. Let 
us introduce one and two particle inclusive multiplicity densities p\ and p 2 based on the inclusive 
differential cross section with respect to the total inelastic cross section o ine \ as follows[]|] 

' do = p\{r])dr] 



Ginel 



d z a = p 2 {rii,ri2)dri l dri2. (3.1) 

With these densities, a two particle density correlation function is defined as 

C 2 (T]I,T?2) =P2(T?l,T?2)-pi0?l)Pl( T ?2)- ( 3 -2) 

The mathematical connection between second order normalized factorial moment F 2 and the two 
particle correlation function is expressed as[^] 

F (8n) = ^ = P2{r)\,r)2)dr)idr}2 

i r r 5r > c 2 (tji,tj 2 ) 



(Sri) 2 JJ P -y 



2 



drjidri2 + l (3.3) 



where n is the number of produced particles, 5t] is the rapidity interval which defines the measur- 
able range of |rji — TJ2I, P\ is the average number density per unit length within Sr\ which is defined 

as 

l r Sr i 

ft = 8jjj Pi(n)dfl- (3-4) 
The two particle correlation function C2 can be parametrized based on the one dimensional 



function form obtained in Eq.( |2.8| ). However, one has to bear in mind that the damping behavior 
in Eq^l^) is originated only from the spatial inhomogeneity of the system in a fixed temperature. 
In many experimental conditions, the initial system temperature can not be specified as a point. 
For instance, corresponding temperature is indirectly discussed by relating it with the collision 
centrality. The centrality bin has a finite size and it causes fluctuations originating from the finite 
temperature bin size. In principle this kind of fluctuations must be independent of the spatial 
fluctuations. In addition, although a self correlation at the zero distance between the two sub 
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volumes in Eq.(2.6) was excluded, the self correlation can not be excluded in the integrated two 
particle correlation function contained in Eq.(fO|), since there is no explicit distinction between 
two particles in the following analysis procedure. Therefore the general function form for the 
normalized two particle correlation function in the one dimensional analysis can be parameterized 
as follows which explicitly contains a constant term j8 ; 

C ^i) =ae -*,/t +p ( 3. 5) 
Pi 

where p\ is proportional to the mean multiplicity in each centrality. 

Instead of F 2 itself, we will use an indirect parameter k of the Negative Binomial Distribu- 
tion(NBD) in the following analysis which is defined as 

Fk ^ n) -T(n-l)T(k) \\+n/k J l+n/k' ^ 

where pL corresponds to the mean value and k~ l reflects deviations from the completely random 
case i.e. the Poisson distribution which corresponds to k = oo. Intuitively k~ l is a measure how 



strongly particles are correlated. The mathematical relation between k and F2 is expressed as[10] 



k~ 1 =F 2 -l. (3.7) 

The reason why we adopt NBD rather than F2 is that NBD can provide an approximate probability 
distribution which enables us to estimate how inefficient or dead areas of the detector system bias 
the k parameter and to obtain the true value of k based on the estimation, while factorial moment 
itself does not provide any specific models on the distribution function which resulted the observed 
factorial moment. 

As the result, the relation between the NBD k parameter and the pseudo-rapidity interval 8rj 



for the parametrization given in Eq.( |3.5| ) is expressed as 

(3.8, 

Once a, £, and /I are obtained from the NBD fits, one can measure the product of the static 



Xk=oT « pi 2 £a = (-0- ) £a (3.9) 



susceptibility and the corresponding temperature based on Eq.( |2.11| ) and (p.5|); 

2 

M_ } 

where \i max is the mean multiplicity in the most central collision event sample in the following 
analysis. Experimentally it is enough to see how the Xk=oT behaves as a function of a quantity 
which reflects T. 

4. Data analysis 

In the data analysis 258k events taken by a minimum bias trigger during Run2 were used. 
The trigger requires coincident hits in Beam-Beam Counters(BBC) and Zero-Degree Calorimeters. 
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Figure 1: Multiplicity distributions in each 8rj indicated inside the figure measured in 0-10% centrality bin 
in Au+Au collisions at ^/snn = 200GeV. The horizontal axis is normalized by the mean multiplicities. The 
vertical axis is scaled by the factors indicated inside the figure. 



The charged particles were reconstructed by a drift chamber and two multi-wire chambers with pad 
readouts without a magnetic field in order to statistically enhance the low p t charged particles. It 
is essential to focus on the low p t charged pions for the discussion of the phase transition. The 
charged particles were measured by the east arm of the PHENIX detector which covers ±0.35 in 
pseudo-rapidity and n/2 in azimuthal angle (j) around the beam line. Events with collision points 
within ±5cm from the center of the detector along the beam line were analyzed. The collision 
centralities were classified by the correlation between BBC charge sum and ZDC energy sum[Q]. 
We have checked the detector stability rigorously for a run range we have analyzed. 

In addition, the dead or inefficient areas in the detector have been identified and the effects 
on the NBD parameters were carefully studied. We have estimated the relation between input k 
values of NBD and biased k values due to dead or inefficient areas for each size of pseudo-rapidity 
interval in each centrality. As long as the baseline distribution is known as NBD, one can generate 
the number of particles based on the probability distribution in all fine i~i — (j) bins respectively. 
Dead maps can be produced from the track projection points in tj — <p plane in the real data and the 
maps were used to suppress particles which hit the dead areas defined in the maps. After applying 
the dead maps, the biased multiplicity distributions were fit with NBD in given rapidity intervals 
again, which provides correction factors between the true k values and the biased ones. Therefore 
systematic errors arise from the definition of the dead maps. As a general tendency, the dead maps 
scales only the absolute values of k, but not change the correlations between k and rapidity interval 
sizes drastically. 
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5. Results 

Fig.[l| shows the charged particle multiplicity distributions in each pseudo-rapidity interval 
from 1/8 to 8/8 of the full rapidity coverage of |tj| < 0.35 with 0-10% events in the collision 
centrality. The distributions are shown as a function of the number of tracks n normalized to the 
mean multiplicity (n). The vertical error bars show the statistical errors. The solid curves were 
determined by performing the NBD fit. In the following analysis, we have performed the NBD 
fit in each pseudo-rapidity interval size from 1/32 to 32/32 of the full rapidity coverage of 0.7 to 
determine a function shape in k vs. 8f] more precisely. The mean and RMS of reduced % 2 values in 
the NBD fit over all centralities and all interval sizes were obtained as 0.69 and 0.72 respectively, 
which corresponds to typically 95% confidence level. Therefore it is good enough to assume NBD 
as a baseline multiplicity distribution to obtain the integrated correlation function through the k 



parameter based on Eq.(3.8). 



Fig.|| shows corrected k parameters as a function of pseudo-rapidity interval sizes for centrality 
classes indicated inside the figure. The upper and lower two panels correspond to 10% and 5% 
centrality bin width cases, respectively. The vertical error bars show the statistical errors and boxes 
show the systematic errors which come from correction factors on k due to the possible variation 
of dead or inefficient areas in the tracking detector. The solid line indicates the fit result by using 



Eq.(p.8|) with errors of quadratic sum of the statistical and systematic errors. The fit was performed 
from 0.02 to 0.7 in pseudo-rapidity. The lowest centrality bin was determined as 55-65%. The fits 
are remarkably well resulting reduced % 2 of 0.44 at the worst which corresponds to above 99% 
confidence level. This guarantees that the parametrization is actually reasonable. 

Fig.[| a), b) and c) show obtained fit parameters a, j3 and t, as a function of the number of 
participants N p where results for both 10% and 5% centrality bin width cases are plotted as red and 
blue circles respectively. N p was obtained from the centrality classes based on the Glauber model 



which is explained in [12] in detail. The horizontal errors correspond to ambiguities on the mean 
values of N p when the centralities are mapped upon N p . The vertical error bars are obtained from 
errors on the fitting parameter by the Minuit program. 

Fig.|| shows the product of static susceptibility and corresponding temperature Xk=oT as a 
function of the number of participants N p in the case of 10% centrality bin width. This quantity is 



proportional to pi at, based on Eq.(2. 11), where p\ is normalized to 1.0 in the 0-10% centrality 



as defined in Eq.(3.9). The horizontal errors correspond to ambiguities on the mean values of N p 
when the centralities are mapped upon N p . The errors on %k=oT were estimated by taking the error 
propagation between the three parameters into account. 

6. Discussions 

Since the a parameter corresponds to the correlation strength at the zero rapidity interval, the 
correlation is expected to be dominated by the Bose-Einstein effect in a short range as demonstrated 



by [ 1 1 ]. Since the observed chaoticity parameter X is almost constant as a function of N p in Au+Au 
collisions at \/Snn =200GeV, the assumption of a constant a over all centrality samples might be 
more reasonable rather than leaving a as a free parameter within the accessible 8rj region. This is 



because Eq.(3. 8) indicates that k can be approximated as l/(2a^/8rj +j3) in the limit of £ > j8r\ « 
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Figure 2: Corrected A; parameters as a function of pseudo-rapidity interval sizes for centrality classes in- 
dicated inside the figure. The upper and lower two panels correspond to 10% and 5% centrality bin width 
cases, respectively. The vertical error bars show the statistical errors and boxes show the systematic errors 
which come from correction factors on k due to the possible variation of dead or inefficient areas in the 
tracking detector. The solid line indicates the fit result by using Eq.(3.8) with errors of quadratic sum of 
the statistical and systematic errors. The fit was performed from 0.02 to 0.7 in pseudo-rapidity. The lowest 
centrality bin was determined as 55-65%. 
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Figure 3: Extracted fit parameters a, j3 and % as a function of the number of participants N p where results 
for both 10% and 5% centrality bin width cases are plotted as red and blue circles respectively. N p was 
obtained from the centrality classes based on the Glauber model which is explained in jl2| ] in detail. The 
horizontal errors correspond to ambiguities on the mean values of N p when the centralities are mapped upon 
N p . The vertical error bars are obtained from errors on the fitting parameter by the Minuit program. 
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Figure 4: The product of static susceptibility and corresponding temperature Xk=oT as a function of the 
number of participants N p in the case of 10% centrality bin width. This quantity is proportional to P\ 2 olE, 
based on Eq.(2.11), where f>\ 2 is normalized to 1.0 in the 0-10% centrality as defined in Eq.(3.9). The 
horizontal errors correspond to ambiguities on the mean values of N p when the centralities are mapped upon 
N p . The errors on Xk=oT were estimated by taking the error propagation between the three parameters into 
account. 



1 and there is an inevitable correlation between a and £, in the limit. Hence, unless the data points 
in very small 5t] region are available, one can not completely neglect the correlation between a 
and in the fit with the integrated two particle correlation function. The direct measurement of 
two particle correlation function as demonstrated by [11] would resolve this issue in future. 

It should be emphasized that the parametrization in Eq.(3.8) is robust and practically works, 
since the j3 parameter can absorb the effect of centrality bin width effects. Fig.|| b) shows the shift 
of N p dependence to lower values from 10% to 5% centrality bin width case, while other physically 
crucial parameters a and ^ are rather stable. In other words, the j3 parameter contains the effect 
from the fluctuations on N p . The ambiguity on N p measured by PHENIX is not large compared 
to e.g. NA49 where only spectators from the projectile nucleus are measurable and it causes an 
increase of scaled variance of multiplicity distributions in peripheral collisions due to dominantly 
large N p fluctuations in the target nucleus. This is due to the partial sampling with respect to the 
total number of nucleons in two colliding nuclei. Since both projectile and target nuclei in both 
sides can be measured by BBC and ZDC at PHENIX, this kind of large ambiguities on N p is more 
suppressed even in peripheral collisions. Even if N p fluctuation is still remaining, j8 parameter 
can absorb this kind of offset parts of fluctuations and the N p fluctuation is not harmful for the 
measurement of correlation lengths at all, since correlation lengths are based on the differential of 
fluctuations. In addition, j8 would contain effects from the azimuthal correlation such as elliptic 
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flows, since the PHENIX does not cover the full azimuth, fluctuations caused by reaction plane 
rotations and elliptic flows should be contained in principle. However, since the flow effect looks 



constant within the experimental accuracy in the pseudo-rapidity direction over unit rapidity[|13|], 
the measured correlation lengths in the pseudo-rapidity direction is not expected to be affected by 
elliptic flows or initial geometrical biases owing to the j8 parameter. 

Concerning the relation between N p and temperature, it is natural to assume that N p can be a 
monotonically and continuously increasing function of the initial temperature at least in the mid- 
rapidity region, since dEj / dr\ per participant pair is slowly raising as a function of N p as already 
observed in several RHIC energies [|]]. The increase of the correlation length is seen at N p ~ 100 
and the corresponding energy density based on the Bjorken picture is eg/T ~ 2.5GeVfm~ 2 which 
has been measured by PHENIX[Q|. It is interesting to note that the energy density coincides with 
the one where the first drop ofj/y suppression from the normal nuclear absorption was observed at 



SPSJ141], though the consistency of Ej scales between two experiments must be carefully checked. 



7. Summary 

The multiplicity distributions measured in Au+Au collisions at ^fs~m = 200GeV are found to 
be well described by the negative binomial distribution. The two point correlation lengths have 
been extracted based on the function form by relating pseudo-rapidity density fluctuations and 
the Ginzburg-Landau theory up to the second order term in the free energy with the scalar order 
parameter. The function form can fit k vs. 8r\ in all centralities remarkably well. The correlation 
lengths as a function of the number of participants N p indicate a non monotonic increase at around 
N p = 100. This could be a symptom of a critical behavior. Within the present systematic errors, 
the product of the static susceptibility Xk=o and the corresponding temperature T does not show an 
obvious turning point at the same N p where the correlation length increases. 
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